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Recent studies have revealed key roles of noncoding RNAs in sex-related pathways, but little is known about the evo- 
lutionary forces acting on these noncoding RNAs. Profiling the transcriptome of Drosophila melanogaster mth whole-genome 
tiling arrays found that 15% of male-biased transcribed fragments are intergenic noncoding RNAs (incRNAs), suggesting 
a potentially important role for incRNAs in sex-related biological processes. Statistical analysis revealed a paucity of male- 
biased incRNAs and coding genes on the X chromosome, suggesting that similar evolutionary forces could be affecting the 
genomic organization of both coding and noncoding genes. Expression profiling across germline and somatic tissues 
further suggested that both male meiotic sex chromosome inactivation (MSG) and sexual antagonism could contribute to 
the chromosomal distribution of male-biased incRNAs. Comparative sequence analysis showed that the evolutionary age 
of male-biased incRNAs is a significant predictor of their chromosomal locations. In addition to identifying abundant sex- 
biased incRNAs in the fly genome, our work unveils a global picture of the complex interplay between noncoding RNAs 
and sexual chromosome evolution. 



[Supplemental material is available for this article.] 

Sex chromosomes are major targets for sex-related selection, and 
many genome-wide studies have revealed differences between sex 
chromosomes and autosomes with respect to divergence rate, gene 
content, and gene expression patterns (Vicoso and Charlesworth 
2006; Ellegren and Parsch 2007; Mank 2009; Qvarnstrom and 
Bailey 2009). Male-biased coding genes, i.e., genes that are more 
highly expressed in males than in females, are unevenly distrib- 
uted between the sex chromosomes and autosomes in Drosophila, 
mammals, and worms (Betran et al. 2002; Parisi et al. 2003; Ranz 
et al. 2003; Khil et al. 2004; Reinke et al. 2004; Wang et al. 2005). In 
mammals, coding genes expressed in male meiotic or post-meiotic 
cells are underrepresented on the X chromosome, whereas coding 
genes expressed in premeiotic stem cells are overrepresented on 
the X chromosome (Khil et al. 2004). Similar results were obtained 
from transcriptional profiling of samples enriched with D. mela- 
nogaster meiotic cells (Vibranovski et al. 2009a). In D. melanogaster 
and in Caenorhabditis elegans gonads, male-biased coding genes are 
found to be underrepresented on the X chromosome (Betran et al. 
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2002; Parisi et al. 2003; Ranz et al. 2003; Reinke et al. 2004). In 
Drosophila this underrepresentation is observed only for old 
X-linked coding genes. Young male-biased coding genes, those 
that emerged after the split of the melanogaster subgroup (<3-6 
million yr ago) (Russo et al. 1995), are found to be overrepresented 
on the X chromosome (Zhang et al. 2010a). 

Directional movement of male-biased genes out of the X 
chromosome is one evolutionary process that could contribute to 
such an uneven chromosomal distribution of male-biased genes in 
these taxa. New Drosophila retrogenes tend to escape from the X 
chromosome and are more likely to be expressed in testis (Betran 
et al. 2002), and excessive male-biased retrogene traffic has been 
observed on the mammalian X chromosome (Emerson et al. 2004). 
Further studies showed that new DNA-based duplicate coding 
genes exhibit a similar chromosomal distribution pattern to ret- 
rogenes (Betran et al. 2002; Emerson et al. 2004; Meisel et al. 2009; 
Vibranovski et al. 2009b), suggesting that the uneven chromo- 
somal distribution of male-biased genes might not depend on 
a specific molecular mechanism but rather is the product of natural 
selection acting on genes with male-related functions (Meisel et al. 
2009; Vibranovski et al. 2009b). This hypothesis is supported by 
independent evidence from population genomic analysis of copy 
number variation of Drosophila retrogenes (Schrider et al. 2011). 
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Several different hypotheses invoking natural selection could 
explain the paucity of X-linked male-biased genes. First, in- 
activation of X-linked genes during male meiosis (meiotic sex 
chromosome inactivation [MSCI]) (Lifschytz and Lindsley 1972) 
may favor the accumulation of testis-expressed genes in autosomes 
(Betran et al. 2002; Emerson et al. 2004). Second, the sexual an- 
tagonism hypothesis predicts that the chromosomal distribution 
of sex-biased genes is driven by sexually antagonistic forces (Rice 
1984), such that dominant alleles with beneficial fitness effects in 
males, but detrimental effects in females, have a higher probability 
of being fixed on the autosomes (Charlesworth et al. 1987). An- 
other dominance-independent hypothesis of sexual antagonism 
driving germline X-inactivation predicted the demasculinization 
of the X chromosome based on different sojourning times of X, 
1/3 of the time in males and 2/3 of the time in females, which 
would result in an excess of male-biased genes out of the X chro- 
mosome and enrichment of X-linked female-biased genes (Wu and 
Xu 2003). Third, the dosage compensation hypothesis predicts 
that hypertranscription of the Drosophila X chromosome in males 
prevents further up-regulation of X-linked male-biased genes, thus 
favoring their relocation to an autosome (Vicoso and Charlesworth 
2009; Bachtrog et al. 2010). Although evidence for all three hy- 
potheses has been demonstrated using coding genes in Drosophila 
(Parisi et al. 2003; Hense et al. 2007; Vibranovski et al. 2009a, 
Bachtrog et al. 2010; Kemkemer et al. 2011, 2013), no systematic 
experimental study on noncoding RNAs has been done so far to 
test these hypotheses. Noncoding RNAs (ncRNAs) play important 
roles in many reproductive processes (Mattick and Makunin 2005; 
Prasanth and Spector 2007). If selection governs the chromosomal 
distribution of sex-biased genes, we expect male-biased ncRNAs to 
exhibit a chromosomal distribution similar to that observed for 
coding genes. 

In this report, we tested these hypotheses by experimentally 
identifying male-biased ncRNAs in D. melanogaster and analyzing 
their chromosomal distribution. Whole-transcriptome profiling 
revealed a large number of intergenic noncoding RNAs (incRNAs) 
with male-biased expression in both whole body and reproductive 
organs, which we confirmed with RT-PCR. We demonstrate that 
these incRNAs are unevenly distributed between the autosomes 
and X chromosome. Comparisons of germline and somatic tissue 
transcriptional profiles suggest that sexual antagonism and male 
germline MSCI both could be contributing to the peculiar chro- 
mosomal distributions of male-biased incRNAs. In concordance 
with previous studies on coding genes, comparative genomics 
analyses revealed that male-biased incRNAs that originated during 
different evolutionary periods have different chromosomal distri- 
bution patterns, indicating that evolutionary time has a significant 
effect on their chromosomal locations (Zhang et al. 2010a). As for 
coding genes (Zhang et al. 2010a), we found that old male-biased 
incRNAs (>6 my old) are enriched on autosomes, whereas new 
male-biased incRNAs are enriched on the X chromosome. 

In addition, our analyses shed some light in the current de- 
bate about the demasculinization of the X chromosome in Dro- 
sophila. More specifically, our analyses clarify why recent studies 
have shown that Drosophila testis-specific genes are not un- 
derrepresented in the X chromosome, but male-biased genes are 
(Meiklejohn and Presgraves 2012; Meisel et al. 2012). Gene age is 
positively correlated with expression breadth (Zhang et al. 2012). 
Therefore testis-specific genes (narrowly expressed genes) are 
enriched with very young genes, which were previously shown to 
be overrepresented in the X chromosome (Zhang et al. 2010a). In 
order to evaluate the chromosomal distribution of testis-specific 



genes in an unbiased way, we analyzed the old testis-specific cod- 
ing genes and found that they are underrepresented in the X 
chromosome. These results corroborate our previous findings that 
the process of desmasculinization is an evolutionary process that 
appears only over evolutionary time in both Drosophila and 
mammals (Zhang et al. 2010a,b). 

Results and Discussion 

Transcriptome profiling reveals abundant incRNAs 
with sex-biased expression 

We first profiled the transcriptomes of whole male and female 
D. melanogaster adults using Affymetrix whole-genome tiling ar- 
rays. The arrays used 3,116,816 25 -nt probe pairs to assay tran- 
scription of 109,088,560 bp of repeat-free euchromatic genome. 
We detected a total of 35,884,625 bp (-32.89%) in male and 
32,921,857 bp (-30.18%) in female as being transcribed. Of those 
nucleotides transcribed, 7,738,215 bp (—21.56%) exhibited male- 
biased expression and 3,649,022 bp (—11.08%) exhibited female- 
biased expression with a greater than twofold increase. In addition 
to whole body samples, we profiled the transcriptome of adult 
reproductive tracts (gonads: testis, ovaries; and accessory glands) 
and nonreproductive tracts (gut and thorax). We identified 
29,188,400 bp (-26.76%), 25,316,156 bp (-23.21%), and 
25,843,450 bp (—23.69%) as being transcribed in testis, ovaries, 
and accessory glands, respectively. Of those nucleotides tran- 
scribed, 10,066,819 bp (-34.49%) in testis and 7,633,727 bp 
(—29.54%) in accessory glands were identified as male-biased with 
a greater than twofold increase over expression in ovaries. As 
expected, a relatively low proportion of male-biased transfrags 
( trans cribed fragments; see Methods for further details) was found 
in gut (7.36%-~2,605,069 bp of 35,390,445 bp) and thorax (9.51%; 
3,043,360 bp of 32,001,727 bp) (detailed transfrag counting and 
proportions could be found in Supplemental Fig. SI and Supple- 
mental Table SI). 

According to FlyBase annotation release 5.46, 10%-20% of 
male-biased and female-biased transfrags belong to intergenic re- 
gions (Fig. 1). Analysis of coding potential with Coding Potential 
Calculator (CPC) (Kong et al. 2007) suggested that >98% of these 
intergenic transfrags are truly noncoding transcripts. We now fo- 
cus on these intergenic noncoding RNAs, or incRNAs. Intronic 
noncoding transfrags were excluded from most of our analyses 
because current annotations might not distinguish true intronic 
noncoding RNAs from novel exons of protein-coding genes or 
nondegraded intronic transcripts (we provide the chromosomal 
distribution of intronic ncRNAs in Supplemental Fig. S2). 

In all three comparisons performed using male/female adult 
transcriptome profiles derived from either whole body or re- 
productive tracts organs (Fig. 1A-I), the fraction of incRNAs is 
significantly higher in sex-biased comparisons than in non-sex- 
biased comparisons (P < 2.2 X 10~ 16 , x 2 test) (Fig. 1A,B versus Fig. 
1C; Fig. 1D,E versus Fig. IF; Fig. 1G,H versus Fig. II), suggesting as 
expected that incRNA expression differences between males and 
females tend to be associated with sex-related biological pro- 
cesses. As expected, comparisons between reproductive organs 
(testis and accessory glands versus ovaries) identified signifi- 
cantly more sex-biased incRNAs than in the nonreproductive 
organs comparisons (32%-35% versus 12%-16%; \ 2 test, P < 
2.2 X 10" 16 ) (Supplemental Fig. S3), suggesting that incRNAs, like 
any other transcription unit, are more prone to be involved with 
sex-related functions in the reproductive organs. Two examples 
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Figure 1. Genomic distribution of sex-biased transfrags. Expression profiling was done with 
Affymetrix whole-genome tiling arrays. Exon/intron/intergenic annotations were retrieved from FlyBase 
(version 5.46). Rows represent comparisons of male and female whole body RNA (A-Q, testis versus 
ovary (D-F), accessory gland versus ovary (G-l), male versus female gut (J-L), and male versus female 
thorax (M-O). Columns represent male-biased (4,0,0,/,/^), female-biased (B,E,H,K,N), and non-sex- 
biased expression (C,F,/,/.,0). 



of sex-biased intergenic noncoding RNAs are shown in Supple- 
mental Figure S4. 

Sex-biased incRNAs are nonrandomly distributed between 
the autosomes and X chromosome 

We tested the hypothesis that selection governs the chromosomal 
distribution of sex-biased genes by comparing the distributions of 
sex-biased coding and noncoding RNAs. If male-biased ncRNAs are 
randomly distributed along the chromosomes, there is probably 
no selection forces acting on noncoding regions, and this is a 
unique property of the coding genes. As predicted by our hy- 
pothesis, an excess of autosomal male-biased incRNAs were iden- 
tified in whole body (21% excess) and testis (21%) (Fig. 2) (odds 
ratios [ORs] = 1.21, P < 0.01, Fisher's exact test, for each compar- 
ison). We applied the Fisher's exact test to assess the uneven 



chromosomal distribution using the esti- 
mated odds ratio (OR) as an intuitive 
measurement. Odds ratio is the ratio be- 
tween male-biased genes (autosomal/ 
X-linked) and non-male-biased genes (au- 
tosomal/X-linked). Thus, an odds ratio >1 
indicates male-biased genes are enriched 
on autosomes, and <1 indicates X enrich- 
ment. (See Methods for more details). A 
similar trend was found for coding trans- 
frags, consistent with previous studies 
(Fig. 2; Parisi et al. 2003; Ranz et al. 2003; 
Emerson et al. 2004; Sturgill et al. 2007). 
In contrast, we found that female-biased 
incRNAs are significantly overrepresented 
on the X chromosome in comparisons of 
ovary versus testis and ovary versus ac- 
cessory gland (OR range = 0.77-0.86, P < 
6.53 X 10" 3 ) (Fig. 2B,C). No significant 
departure from random chromosomal 
distribution, however, was observed for 
female-biased incRNAs derived from whole 
body female versus male comparison 
(OR = 1.01, P = 0.44) (Fig. 2A). 

We verified sex-biased expression of 
incRNAs using tissue-specific RT-PCR. We 
previously detected 528 incRNAs in adult 
flies in a random sample of D. melanogaster 
intergenic regions using the RT-PCR ap- 
proach (Li et al. 2009). Using the same 
primers, we detected about half of these 
incRNAs (261 of 528) in testes, ovaries, 
and heads in Drosophila adults (RT-PCR 
primers are listed in Supplemental Table 
S2). We considered testis-biased incRNAs 
those that were detected in testis but not 
in ovaries. The reverse logic was applied 
to ovary-biased incRNAs. We considered 
non-sex-biased incRNAs those that were 
detected in both testis and ovary and those 
that were detected only in head. The 
chromosomal distribution of testis-biased 
and ovary-biased incRNAs differed signifi- 
cantly from the distribution of non-sex- 
biased incRNAs (P = 0.0287, Fisher's exact 
test) (Supplemental Table S3), with a sig- 
nificant deficiency of X-linked testis-biased incRNAs (OR = 0, 
P = 0.0223 for testis versus non-sex-biased incRNAs) and 
a marginal enrichment of X-linked female-biased incRNAs 
(OR = 1.35, n.s., probably due to small sample size). Thus the 
PCR-based independent test revealed the same robust chro- 
mosomal distribution patterns for Drosophila male-biased incRNAs, 
verifying that the observed underrepresentation of male- 
biased incRNAs in the X chromosome is not a methodological 
artifact. 



Disentangling the contribution of sexual antagonism 
and MSCI to the demasculinization of the X chromosome 

X chromosome demasculinization is the evolutionary process by 
which selective forces drive male-biased genes off the X chromo- 
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Figure 2. X chromosome demasculinization was observed for incRNAs as well as coding transfrags, based on comparison of whole body of males versus 
females (A), testes versus ovaries (B), and accessory glands versus ovaries (C). Comparison of the gut and thorax of males versus females is shown in D and 
E, respectively. Odds ratios >1 indicate enrichment on autosomes, and <1 indicates enrichment on X chromosome. Significant deviations are indicated: 
(***) P< 0.001; (**) P< 0.01, Fisher's exact test. Blue and red bars represent male-biased and female-biased transfrags, respectively. 



some, either relocating them to the autosomes or eliminating 
them from the genome entirely. We investigated the contribution 
of MSCI and sexual antagonism to the observed X chromosome 
demasculinization for incRNAs in Drosophila. It is not trivial to 
separate the effects of sexual antagonism and MSCI, as MSCI only 
occurs in male meiosis, but sexual antagonism may occur in any 
tissue or cell. However, most sex-biased expression is found in 
testes and ovaries, especially in the meiotic phase (Parisi et al. 2003; 
Vibranovski et al. 2009a). MSCI could therefore be assessed 
by analyzing testis-expressed genes with biased expression in 
meiosis but not in mitosis, thus including the effect of inactivation 
of X-linked genes in meiotic cells but ignoring sexual antagonistic 
effects present in mitotic cells (Vibranovski et al. 2009a). However, 
meiotic cells could also be under the effect of sexual antagonistic 
forces preventing the complete separation of those two processes. 
We thus identified male-biased RNAs involved in meiosis as those 
that are testis-biased but not accessory-gland biased, as accessory 
glands only contain mitotic cells. Accessory glands produce pro- 
teins and compounds that comprise seminal fluid and affect the 
reproductive capacity of both sexes (Ravi Ram and Wolfner 2007). 
Accessory gland-biased genes are therefore potential sexually 
antagonistic genes. We observed a statistically significant over- 
representation of strictly testis-biased incRNAs on the autosomes 
(OR = 1.20; P = 0.0069), suggesting that MSCI contributes to the 
desmasculinization of the X chromosome despite the effect of 
accessory gland expressed genes with sexual antagonistic effects. It 
is possible that there are sexually antagonistic genes expressed in 
testis mitotic cells that are not expressed in accessory gland mitotic 
cells. Therefore our data only suggest the role of X-inactivation in 
producing the paucity of X-linked male-biased incRNAs. Never- 
theless, the role of MSCI observed for incRNAs was also observed 
for coding exons (OR = 1.50, P = 1.07 X 10~ 31 ) and is consistent 
with previous observations derived from protein-coding genes 
expressed in meiosis (Vibranovski et al. 2009a). 

Conversely, we assessed the effects of sexual antagonism by 
comparing incRNAs with biased expression in accessory gland 



(somatic) but unbiased in testis (spermatogenesis). Statistical tests 
showed no significant X demasculinization for these incRNAs 
(OR = 1.01, P = 0.464). Although we found no X demasculinization 
of accessory gland-biased incRNAs, it should be noted that acces- 
sory gland-biased coding exons are significantly underrepresented 
on the X chromosome even after removing genes that are also 
testis-biased (OR = 1.41, P = 2.82 X 10" 12 ). Moreover, sex-biased 
genes in other male-specific tissues may be sexually antagonistic; 
accessory gland is not the sole male-specific somatic tissue in 
Drosophila (e.g., male genitalia) (Liu et al. 1996). 

To further investigate the effects of sexual antagonism in 
other somatic tissues, we performed additional transcriptome pro- 
filing on nonreproductive organs: the gut and thorax of male and 
female adults. No significant chromosomal distribution imbalance 
for either male-biased or female-biased coding exons was found 
(Fig. 2D,E). Female-biased incRNAs expressed in the thorax and gut 
do not deviate from the random chromosomal distribution in both 
tissues. Male-biased incRNAs are also randomly distributed on 
chromosomes, except those expressed in the gut, which are sig- 
nificantly overrepresented on the X chromosome (Fig. 2D), the 
opposite pattern expected for demasculinization. One possible 
explanation for the excess of X-linked male-biased genes found 
only in the gut is a higher proportion of young genes, which are 
known to be found in excess in the X chromosome (Zhang et al. 
2010a). Indeed, 10% of male-biased genes expressed in the gut 
originated <3 million yr ago in comparison to 7% of those 
expressed in the thorax (Fisher's exact test, P = 0.066). The gener- 
ally random chromosomal distribution of sex-biased genes in the 
gut and thorax, except for the excess of male-biased incRNAs in the 
X, suggests that the demasculinization of genes is more often as- 
sociated with reproductive organs. 

Moreover, our entire data, which combine reproductive and 
nonreproductive organs, suggest that the X demasculinization 
effects of sexual antagonism are limited to accessory gland-biased 
coding genes. These results support the hypothesis that sexual 
antagonism probably contributes less than MSCI to the non- 
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random chromosomal distribution of male-biased genes in 
Drosophila. 

Our results from nonreproductive organs do not support the 
involvement of dosage compensation in generating a paucity of 
male-biased protein-coding genes observed in the Drosophila X 
chromosome (Vicoso and Charlesworth 2009; Bachtrog et al. 
2010). That is, although the thorax and gut also experience dos- 
age compensation, we observed no paucity of X-linked male- 
biased protein-coding genes in those tissues. 



whereas the opposite pattern is found for older male-biased genes 
(Zhang et al. 2010a,b). Therefore, the process of desmasculiniza- 
tion is an evolutionary process that appears over evolutionary 
time. In order to evaluate the chromosomal distribution of testis- 
specific genes in an unbiased way, we should take into account 
only older genes. Moreover, neither study accounted for the re- 
lationship between gene age and expression breadth. Testis-specific 
genes are genes narrowly expressed in testis. Young genes tend to be 
more narrowly expressed than older genes (Zhang et al. 2012). We 



X chromosome demasculinization 

Two recent studies have questioned both 
the demasculinization of the X chromo- 
some in Drosophila and the contribution 
of MSCI to the phenomenon (Meiklejohn 
and Presgraves 2012; Meisel et al. 2012). 
In both papers, the authors claimed that 
a deficit of male-biased genes on the X 
chromosome is attributable solely to 
lower average expression of genes on the 
X relative to the autosomes in Drosophila 
testes, most likely due to an absence of 
dosage compensation in the germline 
(Meiklejohn and Presgraves 2012; Meisel 
et al. 2012). In both studies, this argu- 
ment is mainly based on the random 
chromosomal distribution of testis-spe- 
cific genes (i.e., those that are expressed 
in testis but not in any other tissue) 
as opposed to the deficit of testis-biased 
genes (i.e., those that are expressed more 
in testis than ovary) on the X chromo- 
some. The interpretations of those results 
are the following: (1) MSCI is not a fac- 
tor that contributes to the nonrandom 
chromosomal distribution testis-biased 
genes because the X chromosome shows 
no deficit of testis-specific genes that, in 
theory, should be under meiotic inac- 
tivation; (2) deficit of testis-biased genes 
should be attributed to lack of dosage 
compensation in the male germline as 
comparisons between testis and ovaries 
are the only analysis that present the 
deficit; and (3) direct comparisons be- 
tween testis and ovary expression do 
not control for correlation between ex- 
pression breadth and sex-biased expres- 
sion. Therefore, by comparing those sex- 
specific tissues, one might obtain evidence 
for the paucity of X-linked testis-expressed 
genes. 

However, those studies did not take 
into account the importance of gene age 
when looking to the random chromo- 
somal distribution of testis-specific genes 
(Meiklejohn and Presgraves 2012; Meisel 
et al. 2012). It is now well established that 
young male-biased genes both in Dro- 
sophila and mammals tend to be more 
frequently found in the X chromosome, 
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Figure 3. Percentage of D. melanogaster X-linked genes that are broadly and narrowly expressed. 
Following the methods in Meisel et al. (201 2), genes narrowly expressed are also called specific genes 
(t > 0.7) in each of four sex-limited tissues, with testis-specific expression and detectable in the sperm 
proteome (testis-SP), narrowly expressed in any of 1 4 tissues (narrow), narrowly expressed in one of 1 0 
non-sex-limited tissues (no sex), or broadly expressed genes (t < 0.4). Percentage of X-linked in the 
genome are shown by dashed lines. Significant deviations (Fisher's exact test) are indicated: (***) P < 
0.001 ; (**) P < 0.01 ; (*) P < 0.02. All X-linked genes (A) were separated according to their evolutionary 
age. New (B) and old (C) genes were defined according to Zhang et al. (201 0a), in which old genes are 
at least as old as the split between Sophophora and Drosophila subgenera. At first (A), there is no paucity 
of testis-specific genes in the X chromosome. However, opposite patterns are found for old and new 
genes: Although old testis-specific genes are underrepresented (C), new testis-specific genes are found 
in excess in the X chromosome (B). 
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tested the hypothesis that the random chromosomal distribution 
of testis-specific, but not of testis-biased, genes is caused by a large 
number of testis-specific genes being newly evolved genes. We 
classified coding genes according to their expression breadth 
following Meisel et al. (2012) and separated older and young 
genes according to ages as defined in Zhang et al. (2010a), 
branch 0 and branches 1-6, respectively. First, we confirmed 
that narrowly expressed coding genes are enriched with new 
genes (Supplemental Fig. S6; Zhang et al. 2012). This pattern is 
also true for testis-specific coding genes (Supplemental Fig. S6). 
Second, young testis-specific coding genes are enriched on the 
X chromosome, whereas older testis-specific coding genes are de- 
ficient from the X chromosome (Fig. 3B,C, respectively). We 
therefore conclude that the result of random distribution of testis- 
specific coding genes (Meiklejohn and Presgraves 2012; Meisel 
et al. 2012) is a consequence of the enrichment of testis-specific 
coding genes with recently evolved coding genes in a short evo- 
lutionary period. Therefore, neither demasculinization nor MSCI 
can be ruled out as important players in determining the chro- 
mosomal distribution of male-biased coding genes in Drosophila as 
older testis-specific coding genes are underrepresented in the X 
chromosome (Fig. 3C). 

Although there is no argument over the presence of 
X chromosome reduced expression in male germline cells, there 
are different opinions for the period of this expression re- 
duction. MSCI studies presented evidence of expression re- 
duction in the meiotic stage (e.g., Vibranovski et al. 2009a), 
whereas one study believed that the reduction is also extended 
to the mitotic stage of spermatogenesis (Meiklejohn et al. 2011). 
However, recent analyses reported three independent lines of 
evidence in favor of MSCI analyzing the expression of testis- 
specific promoter reporters, testis from larval stages, and from 
meiotic arrest mutants (Deng et al. 201 1; Vibranovski et al. 2012; 
Kemkemer et al. 2013). Nevertheless, the work presented here 
does not provide evidence in favor or against MSCI, but the 
patterns found are consistent with the phenomenon. 

Male-biased incRNA gene linkage depends on gene age 

Both male-biased incRNAs and male-biased coding transfrags are 
significantly deficient from the X chromosome, but this trend is 
stronger for coding than for noncoding transfrags (Fig. 2). The 
comparison of testis-biased coding vs. noncoding transfrags shows 
OR = 1.21 versus 1.49 (P= 6.63 X 10" 4 ), and the entire gene unit vs. 
noncoding transfrags shows OR = 1.21 versus 1.39 (P = 0.026) (see 
Supplemental Table S4 for details). Do different evolutionary ages 
of male-biased coding and noncoding genes play any role in 
determining the evolutionary dynamics? We inferred the evo- 
lutionary age of incRNAs through comparative sequence anal- 
ysis of the 12 sequenced Drosophila genomes (Drosophila 12 Ge- 
nomes Consortium 2007; Stark et al. 2007). Given the relatively fast 
evolutionary rate of incRNAs (Pang et al. 2006), we took a conser- 
vative dating strategy (see Methods). Following the parsimony 
principle, 23,165 incRNAs were assigned to a unique phylogenetic 
branch compared to out-group species. Among those, 2660 (or 
11.48% of 23,165) were identified as "male-biased" in D. mela- 
nogaster during at least one comparison (Fig. 4A). By comparing 
the chromosomal distribution of incRNAs across two age 
groups, i.e., before and after the split of the melanogaster sub- 
group (3-6 million yr ago) (Fig. 4), we found that older male- 
biased incRNAs, those that originated before the split, are sig- 
nificantly enriched on autosomes (OR = 1.18, P = 0.037 for 



whole body; and OR = 1.38, P = 3.096 X 10" 4 for gonad). In 
contrast, young male-biased incRNAs show the opposite pat- 
tern in both whole body (OR = 0.72, P = 0.05) and gonad (OR = 
0.71, P = 0.03) comparisons. Furthermore, we found that a 
significantly larger percentage (10.86%, or 289 of 2660) of male- 
biased incRNAs emerged very recently (<3 million yr, on 
branches 5 and 6 in Fig. 4A) compared to male-biased coding 
genes that emerged during the same period (4.03%, or 400 of 
9931) (Fisher's exact test, OR = 2.90, P < 2.2 X 10" 16 ). A sig- 
nificantly negative correlation (Spearman correlation rho = 
0.95, P = 0.0008) between the age of the lineages and the 
proportion of X-linked male-biased incRNAs was observed 
(Fig. 4B). 

The age analyses implemented in this study work on the DNA 
sequence level. That means, for sex-biased and unbiased transfrags, 
we can infer the age of the corresponding DNA locus based on the 
presence or absence of information. However, although the DNA 
sequence can be old, the transcription pattern may have only re- 
cently evolved. In this sense, our strategy provides an upper age 
estimate for the expression pattern. Since male-biased expression 
has a higher turnover rate, such an approximation may be too 
conservative, and the age of male-biased transfrags could have 
been overestimated. Therefore, there could be an even larger pro- 
portion of younger male-biased incRNAs, further strengthening 
our conclusions. 

Zhang et al. (2010a) reported an excess of X-linked new pro- 
tein coding genes in Drosophila that had been recently generated 
from DNA-level duplication or de novo gene origination, and the 
proportion of male-biased genes among the X-linked new genes 
diminishes with gene age. The incRNA genes in this study show 
a similar pattern to the new protein-coding genes. However, it 
appears that turnover is more recent in the incRNA genes, con- 
sistent with the known rapid evolution of ncRNA genes (Pang et al. 
2006) and higher turnover rate of microRNA genes (Lu et al. 2008). 
These data indicate that different evolutionary forces, e.g., MSCI 
and sexual antagonism, might play roles at different evolutionary 
timescales (Zhang et al. 2010a). 

In summary, we experimentally identified male-biased non- 
coding RNAs in D. melanogaster and analyzed their chromosomal 
distribution. The identification of a large number of incRNAs that 
showed male-biased expression patterns may explain the signals of 
natural selection previously detected in the noncoding genomic 
regions (Andolfatto 2005). By systematically profiling the whole 
transcriptome of D. melanogaster male and female adult whole 
bodies as well as reproductive tract organs, we revealed a long-term 
removal of male-biased incRNA genes from the X chromosome 
resulting in an uneven distribution of male-biased incRNAs between 
X and autosomes. This led to the long-term X chromosome de- 
masculinization, probably through sexual antagonism and MSCI. 
Finally, we identified distinctive chromosomal preferences between 
young and old male-biased incRNAs. This pattern of male-biased 
incRNAs further generalizes the uneven pattern of male-biased gene 
content on Drosophila autosomes and X chromosomes and suggests 
that the pattern is shaped by natural selection acting on male 
functions. Our results contribute to a global picture of sex chro- 
mosome evolution in the Drosophila genome. 

Methods 

Data sources 

Genome sequence and annotations were obtained from the UCSC 
Genome Browser database (http://genome.ucsc.edu; dm2). Gene 
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model annotations were obtained from FlyBase (http://www. 
flybase.org; v5.46, downloaded in March, 2013). The Affymetrix 
array BPMAP annotation was downloaded from the Affymetrix 
website (http://www.affymetrix.com/products_services/arrays/ 
specific/ drosophila_tilingl_Or.affx). The 12 sequenced Drosophila 
genome sequences (CAF1) were downloaded from http://rana. 
lbl.gov/drosophila/. 

Sample preparation and microarray hybridization 

We extracted total RNA from whole bodies, thoraxes, digestive 
tracts (guts), testes, accessory glands, and ovaries from virgin 



Oregon R adults using the Qiagen RNeasy Mini Kit with on- 
column DNase digestion. All flies were virgin and aged for 1-6 
d post-eclosion before being used in extractions or dissections. 
Whole body and thorax RNA was immediately extracted from 50 
males or females. Gut tissue was completely removed from tho- 
raxes. Testes, accessory glands, ovaries, and guts were dissected, 
placed in RN Alater (Qiagen), and stored at -20°C until RNA ex- 
traction. Roughly 200 testes, accessory glands or ovaries, or 100 
guts were required for each replicate. Testes were separated from all 
other reproductive tract tissues (seminal vesicles, accessory glands, 
and ducts). For statistical independence, all tissues were harvested 
from independent sets of flies. Each of the biological replicates 
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was labeled with GeneChip WT Sense Target Labeling and Con- 
trol Reagents (Part#: 900652) and hybridized to an Affymetrix 
D. melanogaster genome tiling array as previously reported by 
Manak et al. (2006). Labeling and hybridizations were performed 
at the University of Chicago Functional Genomics Facility. Three 
replicates of hybridization were performed for each tissue. 

Tiling array data processing and analysis 

Affymetrix Tiling Analysis Software (TAS vl.1.02) was used to 
process raw tiling array data (Cheng et al. 2005; Manak et al. 2006). 
Raw data were normalized by quantile normalization and the 
median of target intensities was scaled to 100. As suggested in the 
Affymetrix TAS user manual (http://www.affymetrix.com/support/ 
developer/downloads/TilingArrayTools/index.affx), each probe po- 
sition was analyzed in a local smoothing window with bandwidth 
(BW) equal to 50 bp (resulting in a window width of 101 bp) for 
better statistical power. To assess the performance of replicates, 
standard Pearson correlation coefficients between replicates were 
calculated pairwise. The significant correlation (Spearman correla- 
tion rho > 0.98, P-value = 1 X 10~ 5 ) indicated reasonable consis- 
tency between replicated samples. 

A one-tail Wilcoxon signed-rank test on all the probes in the 
window was performed with the alternative hypothesis that the 
true intensity difference between the perfect match (PM) probe 
and mismatched probe is significantly greater than zero. Only 
probes with a P-value < 0.1 were called "positive." Neighboring 
positive probes with max-gap 50 bp and a minimum run of 90 bp 
were grouped as transfrags (transcribed /raiments), then Unified 
Transfrags, or UTS, were further derived by assembling overlapped 
transfrags ("supporting transfrags") in different samples as suggested 
by previous literature (Cheng et al. 2005; Manak et al. 2006). 
"Present" call was produced for each unified transfrag of UTS. A 
unified transfrag will be called "Present" in a given sample if and 
only if at least one of its supporting transfrags is identified in this 
sample. Moreover, the median intensity value of all comprising 
probe sets within each unified transfrag was calculated for each 
sample. Among UTS, —48% unified transfrags are tissue-specifically 
transcribed, and 78% consist of multiple supporting transfrags 
with >50% overlap. Further inspection indicated that 73.52% 
(68,310 of 92,916) of coding exons on autosomes (chr2L, 2R, 3L, 
and 3R) and X chromosome annotated in FlyBase were detected as 
expressed in at least one of the five D. melanogaster samples (adult 
male/female whole body and three reproductive organ tissues 
[testis, ovary and accessory gland]) with >70% coverage. 

An integrated procedure was applied to assess sample-biased 
expression for each unified transfrag (Suppplemental Fig. S5). To be 
conservative, a global Kruskal-Wallis test ("nonparametric one-way 
ANOVA") among the five samples was applied before the pairwise 
Mann-Whitney [/-test. To assess the detection power of our pro- 
cedure, we compared the identified sex-biased protein-coding trans- 
frags to the Sebida database that integrates data derived from multiple 
previous high-throughput studies comparing male versus female 
protein-coding genes expression in D. melanogaster (Gnad and Parsch 
2006). Up to 88% (613 of 695) male-biased and 84% (663 of 789) 
female-biased genes in the Sebida twofold high quality data set were 
also identified by our procedure with > 70% length coverage, sug- 
gesting a high consistency of our procedure with previous studies. 

Annotation of detected transfrags 

Transfrags were classified as "coding" or "intergenic" based on 
their genomic coordinates. Because our work focused on non- 
coding genes, as a conservative estimation, we only consider 
a transfrag as "intergenic" if it does not overlap with any annotated 
FlyBase protein-coding gene models on both the sense and anti- 



sense strands. We further assessed the coding potentials of these 
intergenic transfrags by a SVM-based classifier (Kong et al. 2007); 
the results indicated that >98% intergenic transfrags are truly 
noncoding. After excluding 945 transfrags showing putative cod- 
ing potential at either sense or antisense strand, we classified the 
remaining transfrags as "intergenic noncoding transfrags." The 
intronic transfrags were excluded from follow-up analyses as cur- 
rent annotations might not distinguish true intronic noncoding 
RNAs from novel exons of protein-coding genes or nondegraded 
intronic transcripts. Moreover, in case of potential bias resulting 
from the relatively larger exon number in protein-coding genes 
(when comparing to noncoding RNAs), we re-ran analysis in gene 
level by assigning transfrags to annotated FlyBase genes according 
to the coordination. 

Assessing the relationship between expression breadth 
and gene age for protein-coding genes 

According to Meisel et al. (2012), microarray signal intensities from 
14 adult D. melanogaster tissues were obtained from FlyAtlas 
(Chintapalli et al. 2007). Expression breadth was calculated 
according to the tissue specificity index, t (Yanai et al. 2005). Genes 
were considered as narrowly (tissue-specific) and broadly expressed 
depending on their t value (t > 0. 7 and t < 0.4, respectively). Testis- 
specific genes were considered to be encoded proteins in the sperm 
proteome if they were found in at least one of the two sperm 
proteomes (Dorus et al. 2006; Wasbrough et al. 2010). Gene age 
was obtained by crosslinking CG identifiers with information 
available from Zhang et al. (2010a). 

Comparative genomics analysis to infer evolutionary 
ages along the Drosophila phylogeny 

In silico comparative sequence analysis was performed with all 12 
sequenced Drosophila genomes similar to the procedure reported 
by Sturgill et al. (2007). We ran NCBI BLAST against genomic DNA 
of each species. To handle the relatively low sequence conservation 
of noncoding genes, optimized BLAST parameters were employed 
as suggested in previous literatures (Korf et al. 2003; Freyhult et al. 
2007). A D. melanogaster incRNA is called "absent" in another 
species if there are no hits with E- value < 1 X 10~ 4 and coverage 
>80% found in that species. After making the "present"/"absent" 
call for each incRNA, we dated their origination along the Dro- 
sophila genus phylogenetic tree (Supplemental Fig. S7) following 
the parsimony principle described below. 

IncRNA x is assigned to branch X if and only if it is called 
"present" in all in-group species of branch X and "absent" in all 
out-group species of X. For example, branch 0 includes incRNAs 
that are "present" in all 12 sequenced Drosophila genomes; branch 
4 includes incRNAs that are "present" in D.mel, D.sim, D.sec, D.yac, 



Table 1. Test for chromosomal distribution of transfrags 
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Thus, odds ratio >1 indicates that male-biased genes are enriched on 
autosomes, and <1 indicates X enrichment. All statistical computations 

were performed by R (http://www.r-project.org). Odds Ratio = 



636 Genome Research 

www.genome.org 



Demasculinization of X-Iinked incRNAs in Drosophila 



and D.ere, but "absent" in D.ana, D.pse, D.per, D.wil, D.moj, D.vir, 
and D.gri. 

Assessment of the uneven distribution of the transfrags 
among X chromosome and autosomes 

We applied the Fisher's exact test to assess the uneven chromo- 
somal distribution using the estimated odds ratio (OR) as an in- 
tuitive measurement (see Table 1). 

Independent RT-PCR assay 
RNA extraction 

Adults were collected within 10 d after eclosion. Tissues like 
heads, ovaries, and testes were separately dissected into tubes 
with TRIzol (Invitrogen). RNAs were extracted following TRIzol 
reagent instructions. About 10 jxg sample RNA was mixed with 
20 ijlL RQ1 RNase-free DNase (lunit/jjuL; Promega), 10 uX 10X 
DNase buffer, 2 uX RNase inhibitor HPR1 (Takara) and DEPC 
water (up to 100 uX), and incubated for 3 h at 37°C. Contami- 
nation of RNAs with DNA was ruled out by PCR amplification of 
two pairs of primers for Gapdh2 and 11171a, taking the extracted 
RNAs as template. 

Reverse transcription 

Five tenths micrograms 6-mer random primer (Takara Company) 
and 5 uX dNTPs (2.5 mM, Takara Company) per microgram of RNA 
sample were mixed in a total volume of <15 jjlL in a tube; the tube 
was heated for 5 min to 70°C to melt secondary structure; the tube 
was cooled immediately on ice for at least 3 min to prevent a sec- 
ondary structure from reforming; 5 jjlL M-MLV 5 X Reaction Buffer 
(Promega Company), 5 uX dNTPs (2.5 mM, Takara Company), 0.5 
uX RNase Inhibitor (HPR I, Takara Company), and 1 jjlL M-MLV 
RTase (Promega Company) were added to the annealed primer/ 
template. Then DEPC water was added up to the volume of 25 uX. 
PCR of primer III 71a was conducted to guarantee no genomic DNA 
contamination, and PCR of Gapdh2 was conducted to guarantee the 
quality of cDNA. 

Data access 

The data discussed in this publication have been submitted to the 
NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm. 
nih.gov/geo/) (Edgar et al. 2002) under accession number 
GSE53421. 
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